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Abstract 

A direct numerical algorithm for solving the time-nonlocal non-Markovian master equation in 
the second Born approximation is introduced and the range of utility of this approximation, and 
of the Markov approximation, is analyzed for the traditional dimer system that models excitation 
energy transfer in photosynthesis. Specifically, the coupled integro-differential equations for the 
reduced density matrix are solved by an efficient auxiliary function method in both the energy and 
site representations. In addition to giving exact results to this order, the approach allows us to 
computationally assess the range of the reorganization energy and decay rates of the phonon auto- 
correlation function for which the Markovian Redfield theory and the second order approximation 
is valid. For example, the use of Redfield theory for A > 10 cm -1 in systems like Fenna-Mathews- 
Olson (FMO) type systems is shown to be in error. In addition, analytic inequalities are obtained 
for the regime of validity of the Markov approximation in cases of weak and strong resonance 
coupling, allowing for a quick determination of the utility of the Markovian dynamics in parameter 
regions. Finally, results for the evolution of states in a dimer system, with and without initial 
coherence, are compared in order to assess the role of initial coherences. 
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I. INTRODUCTION 



The quantum mechanics of open systems, i.e. systems interacting with an external en- 
vironment, is currently the focus of widespread attention [l- 3)]. Of particular recent interest 
is the issue of the extent to which quantum coherence of the system is maintained in the 
presence of the environment. Two significant modern examples may be noted: (a) the need 
to maintain coherence in order to implement methods for quantum mechanically controlling 
molecular processes jj], and (b) issues of the role of such quantum coherent processes in 
natural environments, such as the observed long-lived coherent electronic energy transfer 
(EET) in photosynthesis . 

In either of these cases, and in many other examples as well, dynamical evolution of the 
open system provide a significant computational challenge. As such, a variety of approxima- 
tions are often invoked to propagate the system, such as the second Born approximation to 
master equations and the Markov approximation, both of which are the focus of this paper. 
In particular, in this paper we introduce a simple method to solve the second Born quantum 
master equation without doing the Markov approximation on the slowly decaying envelope 
of the density matrix. In addition to being straightforward, this approach allows, by com- 
paring to results using the Markov approximation, a reliable determination of the range of 
coupling strengths and decay rates of the bath auto-correlation function, over which one can 
use the Markovian theory and the second order approximation. In addition, by examining 
the size of the fourth order term, this approach affords an estimate of the range of validity 
of the widely used second-order approximation for model photosynthetic systems. 

Although the method developed here is applicable to general systems with exponential 
bath correlation functions, for computational simplicity we study, as do others, the dimer 
system, generally regarded as a simple photosynthesis EET model. In this case our Marko- 
vian analysis contrasts with, for example, that in Ref. [6| in which Markovian Redfield 
theory is used apriori and its consequences analyzed, as opposed to comparison with exact 
results. 

The particular challenge arising from these EET types of systems relates to the param- 
eter range in which the dynamics occurs. Specifically, quantum dynamics can be readily 
analyzed in two limiting cases defined by the relative contributions of the inter- system cou- 
pling V responsible for excitation energy transfer (EET), and system-bath coupling constant 
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asBi responsible for decoherence. These parameters define two important time scales: the 
excitation transfer time scale Transfer = h/V, and the decoherence time scale Td eC o = Ti/o-sb- 
If the system-bath coupling is very weak and r deco 3> r trans f er , the system is almost closed 
and the Schrodinger equation can be used to study the dynamics. In the opposite case 
Tdeco *C Transfer (strong system-bath coupling), the system is open, the decoherence rate 
is very fast, the dynamics is almost incoherent and a simple Pauli type master equation 
description suffices. These limiting regimes are well understood. Many real systems, such 
as a number of harvesting systems [7Hl8||. however, fall between these extremes. Recent 
observations [5] of the long-lived EET has reactivated interest in these systems. 

The standard approach used to treat this intermediate regime is to use the second Born 
quantum master equation 1|, a perturbative master equation up to second order in system- 
bath interaction with weak system-bath coupling, plus its Markovian approximation (e.g., as 
in the Redfield master equation). Recently, two approaches have been studied for arbitrary 
coupling regimes. One is based on weakening the system-bath coupling removal of system- 
bath interaction and repartitioning the Hamiltonian term using a polaron transformation, 
followed by the standard second Born master equation [3]. The second approach is based 
on a reduced hierarchy equation of Kubo and Tanimura, starting from the path integral 
approach for quantum dissipative systems [20]. Additional methods are also being developed 
by the community. Here, as noted above, we introduce and utilize a particularly direct 
approach. 

In the following section (Section II) we outline the basic model for a dimer. In Section III 
we introduce the second Born quantum master equation and phonon correlation function, 
diagonalize the Hamiltonian and cast the master equation into both the site and energy 
representations. Section IV gives a new auxiliary function method of solving these equations, 
and provides an analysis of results in the Markovian approximation. Discussion of the results 
and the underlying physical picture is given at the end of Section V. In Section VI we discuss 
the regime of validity of the second order approximation in the master equation by estimating 
the order of magnitude of the fourth order term. 

The vast majority of treatments in the literature utilize coherence-free initial conditions 
and study the subsequent dynamics. Results for these initial conditions are compared to 
that obtained with model excitation with weak light in Section VII. The last section provides 
a brief conclusion. 
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II. THE MODEL: DIMER SYSTEM 



Consider a model dimer system given by the following standard Frenkel exciton Hamil- 



tonian jg]: 



H tot = H el + H reora + H ph + H el ~ ph (1) 

2 

H el = ^6»(n| + J(|l)(2| + |2)(l|) (2) 



n=l 
2 



H reor9 = J2 X n\n)(n\, A n = ^2hUi<Pj2 (3) 

n=l i 
2 

ffP h = £ hP h } h ph = £ hUi{p 2 + g 2 )/2 (4) 

n=l i 
2 

H d-ph = ^2v n u n , V n = \n)(n\, u n = -^2huid ni qi (5) 

n=l i 

Here \n) represents the state in which only the n th site is excited and all others are in 
the ground state. The quantity e° is the excited electronic energy of the n th site in the 
absence of phonons, and J is the electronic coupling between the sites which is responsible 
for EET. The ground state energies of the donor and acceptor are set equal to zero and Xj 
is the reorganization energy of the j th site that is dissipated in the bath after the electronic 
transition occurs. The quantity dji is the dimensionless displacement of the equilibrium 
configuration of the i th phonon mode between the ground and the excited electronic state 
of the j th site, and qi,Pi are the dimensionless coordinates and momenta of the i th phonon 
mode of frequency Wj. 



III. THE SECOND-BORN QUANTUM MASTER EQUATION 

The method of projection operators used to obtain open system master equations is well 
known With the help of projection operators one can obtain the following quantum 
master equation for the reduced density matrix of the system in the second Born approxi- 
mation, which is valid when system-bath coupling is weak as compared to the characteristic 
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energy scale of the system [see, e.g., Ref. 



.a. 



dt h 

J 



(Cy(f - r)[y/(t), VjWOO] - C^t - r)[V/(t), /(r)^(r)]) (6) 



Here, the interaction representation has been used, which is defined for system operators as, 

O'it) = Ul(t)OUs(t) (7) 

where Us(t) = exp(—^H s t) is the time evolution operator, and H s = Yln=i( e n + ^n)\n)(n\ + 
J(|l)(2| + |2)(1|) is the system Hamiltonian. Here, the bath is assumed to be a continuum 
of harmonic oscillators, and the bath correlation functions are defined as 

C^^iu^uM)-^)^) ■ (8) 

Below, the canonical average of the bath operators, (uj), which involve the averaging over the 
product of displacement and bath position co-ordinates is taken to be zero. The above master 
equation [Eq. (6)] is also termed the time convolution equation and can be obtained from 
the Nakajima-Zwanzig equation with a zeroth order approximation to the time evolution 
operator in the kernel 

Converting this master equation [Eq. (6)] back to the Schrodinger representation gives 



dp® 



(9) 



dt 

{C l3 (t - rWiit), U s (t - r)V jP (r)Ul(t - r)} - C* 3 (t - r)[V h U s (t - r)p(r)V 3 Ul(t - r)]) . 

We consider the case where the characteristics of the bath as seen by both the sites are 
the same, and there is no bath correlation between the sites. The bath correlation function 
is then of the form Cij(t) = C(t)5ij, where 

/+°° A, , 
x ^ C{U)e ' lUjt - (10) 

C(u) = 2h(l + n(u))J(u), (11) 

where n(ou) is the Bose-Einstein distribution function. Assuming the Drude-Lorentz model 
for the spectral density J(u>) = 2A aj ^ 2 where A is the reorganization energy, and assuming 
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the high temperature approximation (j^pf << 1), as is appropriate for the systems like the 



FMO model [see Refs. 



6) and 20)], we obtain the correlation function as, 
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A. Explicit Site Representation of the Non-Markovian Master Equation 

For the explicit site representation we need the eigensystem of the Hamiltonian. Assuming 
Ai = A 2 = A the eigenvalues Ei and eigenvectors |e$) for the system Hamiltonian 

2 

H s = $>° + X n )\n)(n\ + J(|l)(2| + |2)(1|) (13) 



n=l 



can be easily obtained as 

Ex,2 = + 4 + 2A T + 4 + 2A) 2 - 4(e?e° - J 2 + A(e? + 4) + A 2 )) 

<, 1 , ! =i(A T Vtf+4J=) 1 A = e°-e°. (14) 

Here the column vectors denote components in the site basis, the eigenkets are normalized 
and, since ot\Oi 2 = — 1, they are orthogonal. With lengthy but straightforward calculations, 
Eq. ([9]) for the reduced density operator can be written explicitly in the site representation, 
using Eqs. ( ITOl) and ( ITTT) . and as a set of coupled integro-differential delay equations, 



dj/i(t) A 4 A _ 7 , y* 7r 



dt ft w /3ft 2 jo 

[771 cos(£i 2 (t - r))i/i(r) + m sm(E 12 (t - t))t/ 2 (t)] 

^ - \(t) - J -{1 - 2x{t)) - i^e" 7 ' f dre* 



dt h ,yy 'h y wy /3r 

[-772 sm(E 12 (t - r))y 1 (r) + 773 cos(£ 12 (t - r))y 2 (r) + 2fiy 2 (r)] (15) 



with 77! = 1, r7 2 = - 7 =A_, 773 = ^4^, E 12 = (E, - E 2 )/h = -^p^, fl = 
Here z(t) = p n {t) = (l|p(t)|l> (site), yi (t) = Re[p 12 (t)], and y 2 (t) = Im[p 12 (t)], 
with subscripts denoting the sites. 
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B. Energy Representation of the Non-Markovian Master Equation 



The kets le^a) in Eq. (JT3J) are the eigenstates of the Hamiltonian H s . The equation for 
a general element of the reduced density matrix in energy representation 



P e ab (t) = (e a \p(t)\e b ), 



is obtained from Eqs. (jUJ) and (jTUJ) as 



(16) 



dp e ab {t) 
dt 



i r l 

-iUabP ab - T2 Y] / dr C(t - t) 
71 i,c,d=l J ° 

-C*(t - r) [VrVfe-^^pUr) ~ Vf^e^^Vjj)] 



with u ab = (E a — E b )/h and 



yac 



Results in the energy representation, using Eq. (fl2l) are discussed below. 



(17) 



(18) 



IV. METHOD OF SOLUTION: NON-MARKOVIAN 



To obtain a solution for the non-Markovian case, we first convert the coupled integro- 
differential equations in the site representation [Eq. f JT5|) ] to a larger number of coupled 
ordinary differential equations, a transformation made possible by the exponential form of 
the correlation function. The resultant coupled ordinary differential equations can be numer- 
ically solved easily. This transformation is performed as follows. First, for computational 
simplicity we put r' = 7T in Eq. ( I15p and then yt = t' in the resulting equations, and define 
three auxiliary functions fi(t'): 



hit') = / dr'e 
Jo 

/ a (0 ee f e T 'y 2 (r')dr' 
Jo 



COs[^(t'-T')]£l(T') 

7 



^[^{t* - t')]Ut') 
7 



f 3 (t') ee / dr'e 
'0 



.E 



-r}2 sm[ 



12 



r E 



(f - r')]yi(T') + T} Z cos[— (t' - T')]y 2 {r') 

7 7 



(19) 
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Here, yi(t') = yi(t'/j), jjiit') = 2/2(^/7) an d we a l so define x(t') = x(t' /j). We then obtain 
six coupled ordinary differential equations, three from Eq. (1151) and three from differentiating 
the three auxiliary functions, giving: 



2J 

J'~ fjJ\ „t'~ (4l\ , El2 „t' ~ U'\ I El2 



hit') - e* y x {t') = e* + -H e * r^ 2 (t') - ) hit') 

h(t') = e* W (f), 



/s(f ) - e'vM') = JrmW ~ — V2ie t 'y 1 (t') - (^) ' / 3 (t') , (20) 

7 V 7 / 

where overdots denote derivatives with respect to t'. These equations can be efficiently 
solved numerically. 

For comparison with other studies, results are given below for the particular initial con- 
ditions: Pll (0) = 5(0) = 1, yi(0) = &(0) = 0, ^(0) = / 2 (0) = / 3 (0) = /i(0)/ 3 (0) = 0. 
These initial conditions (corresponding to all the population being on site 1, and no coher- 
ences), are those which have been used extensively in previous investigations [see Ref. jo]] 
but are somewhat unphysical, because they lack initial coherences which become important 
in photo-excitation. We treat this problem of initial conditions and state preparation with 
a more plausible model in Section VII. 



V. ENERGY REPRESENTATION AND MARKOVIAN LIMIT 



A. Formalism 



To consider the Markov approximation, we note that it is particularly simple to invoke 
in the energy representation. Hence, below we first utilize the energy basis and then convert 
the result back to the site representation for comparison with the non-Markovian solution. 

The Markov approximation can be performed when the time scale on which the envelope 
of the density matrix decays is much longer than the decay time of the phonon correlation 
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function [1]. Then one can introduce the following approximation: 

fa(t - r) = e^MjfJt - t) ~ e-^^p^t) = e^p e ab (t). (21) 

As discussed in Ref. [6j, the non-Markovian regime is marked by slow dissipation of the 
reorganization energy (i.e., the slow decay of the phonon correlation function as compared 
to relaxation dynamics time scale, the decay of the envelope part of the density matrix). 
Transitions occur in accord with the vertical Franck-Condon principle. In the Markovian 
regime phonon relaxation is very fast (e.g., large 7) as compared to the decay of the envelope 
of the density matrix. Thus, phonons remain effectively in equilibrium during the EET 
process in the Markovian regime Q|. 

To obtain the equations in the Markov approximation, Eq. (17) is first converted to 
dimensionless form with r' = 7T and if = yt. Putting t — t = r 1 in the resulting equation in 
the energy representation and then implementing the above approximation on the density 
matrix elements allows the time integration to be performed easily for the case of exponential 
phonon correlation function [Eq. (11)]. The result is the set of Markovian equations; 

2X E(^-(0-^«.(o). (-) 

i,c,d 



Here p e ab {t' /^) = p^it'), u ab = Equation fl22|) constitutes a system of coupled ordinary 
differential equations that can be solved with given initial conditions. 

The results can then be transformed back to the site representation using the transfor- 
mation 

Pa{t) = (Mt)\j) = 52<t\ea)fa(eb\j). (23) 

a.b 

where p^ is in site representation, p e ab is in energy representation, and i,j,a,b G {1,2}. 
Equation ( |23l) constitutes four linear equations that provides the relationship between the 
representations. 
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B. Limiting Cases: Analytical Results 

1. Strong Coupling Case: J S> A 

For J ^> A, we have [from Eq. f[T4"|) ] c\\ ~ 1, ct 2 — —1, and ~ 1/2 for i = j and 
~ —1/2 for i ^ j = 1,2) and V^ 3 — 1/2 for all One can then analytically solve 

Eqs. (122]) to obtain the simple expression 

I 4A t i 

Pn(t) = 7>( e f3(4f2+n2 ^ + 1), (24) 

for the traditional initial conditions p\i{t' = 0) = 1, Pi 2 (t' = 0) = p\]i$ = 0) = 0. The 
Markov approximation can be performed when the time scale on which the envelope of the 
density matrix decays is much longer than the decay time of the phonon auto-correlation 
function. Hence, ~pujT^^T\ ^ 1 must hold for the Markov approximation to be valid in 
the J ^> A domain. (Note that the decay time constant for the phonon auto-correlation 
function is unity, since we defined t' = jt.) 

2. Weak Coupling Case: J <C A 

For J < A, we have [from Eq. (TTIj)]. «i, 2 = ^7 (A =f VA 2 + 4 J 2 ) ~ ^(1^1). Hence, 
in this domain a\ ~ 0, and a 2 ^ A/ J. This leads to V^ 1 = V^ 2 = V 21 ~ 0, V 22 ~ 1. 
V 2 n ~ 1, 1/ 2 12 = V 21 ~ J/A, and V 2 22 = (f ) 2 . From Eq. (E2D we then have 



&(0 

2/ A 
ft 2 /3 7 2 



(2(j/A) 2 (r + r*) - 4(j/A) 2 (r + r*)^(0 + 2(J/A)(fo(f) + ^ (*'))) - (25) 



4A 

((J/A)r(2^(0 - 1) - (1 + 2r(J/A) 2 )pl 2 (t') + 2(J/A) 2 r*p* 1 (t')) . (26) 



ft 2 /3 7 2 



where T 



The reduction from a large number (33) of terms in Eq. (|22|) to the smaller number 
of terms [in Eq. fl26]) ] is made possible by neglecting small terms of the order of (J/A) 3 . 
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However, even in this approximation the above equations do not admit a simple analytic 
solution, and no simple analytic expression can be given for the range of validity of the 
Markov approximation. Hence, we invoke a further approximation, neglecting second order 
terms (J/ A) 2 as compared to first order J/ A. By separating real and imaginary parts as 
pf 2 (t') = x{t') + iy(t') and writing p e n (t') = r(t'), Eqs. (ESJ) and (T2§D become 



*v = ^ (f,) - wk y(t,) ~ (27) 

These coupled ordinary differential equations have the straightforward solution: 



[brj - a£]e£cos(77t')e~^' + H + sm(r]t')e~^'), 
x(t') = (acos(t]t') — b sm(rjt')) , 

y {t') = e~^'(asin(r]t') + bcos(rjt')). (28) 

with initial conditions r(t' = 0) = 1, x(t' = 0) = a, y(t' = 0) = b. Here £ = ^r^, r) — j-, 
and e = ^. Interestingly, for a = 6 = the density matrix elements do not change with 
time. This is due to our approximations of only retaining terms first order in J/ A. Thus, for 
the condition of Markov approximation to hold requires £ = "C 1 (again noting that 

the decay time constant for the phonon auto-correlation function is unity since t' = '-ft). 

These analytic results are summarized in Table I, and these inequalities have been nu- 
merically verified (e.g., see Fig. [TJ). Note that the J ^> A result goes over to the J < A 
result as J gets smaller. 



C. Computational Results 

In other parameter regimes, the validity of the Markov approximation [Eq. ( 1221) ] must 
be determined by numerical comparisons with the exact result [Eq. (I20I) ]. Figure [2] com- 
pares the solution for the Markovian master equation to the non-Markovian results for the 
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TABLE I: Regimes of validity of the Markov approximation 



Case 


Approx. matrix elements 


Markovian approximation 


J » A 


i,j,k = 1,2 


4A „ „ i 

/3(4J 2 +ft 2 7 2 ) ^ 1 


J « A 




+ ^«1 



J = 100 cm -1 » A = 1 cm -1 



J = 100 cm -1 » A = 1 cm -1 



( /3 (4 J' 



4 A 



2 + hbar A 2 y A 2) 



« 1 



f " 

V P (4 / A 2 + hb 



bar A 2 y A 2) 



Pll 
1.0 



0.8 
0.6 
0.4 

0.2 1 




10 



15 



20 




15 



20 



7 = 1 cm -1 << A = 100 cm -1 J = 1 cm" 1 « A = 100 cm" 1 



4 A 



hbar A 2 y A 2 



« 1 



4 A 



P hbar A 2 y A 2 



Pll 
1.0000 

0.9999 

0.9998 

0.9997 

0.9996 

0.9995 




20 40 60 80 100 




FIG. 1: Verification of the analytic inequalities (Table I). Time evolution of population on site 
1: blue (solid) curve is the non-Mar kovian solution and red (dotted) curve is the Markovian 
approximation. Upper row (case J >> A): juj^fii^ = 0-04 << 1, A = 2 (upper left graph), 



and 



1A 



/3(4J 2 +?i 2 7 2 ) 



7 



0.9 ~ 1, A = 50(upper right graph). Lower row (case J « A): iA 
0.02 << 1, A = 2 (lower left graph), and A A 2 ~ 3, A = 10 (lower right graph). Clearly, graphs 



are in accord with Table I. Time t' is the dimensionless time, 10 units on this scale are equivalent 
to one ps. 

standard electronic coupling parameter values in photosynthetic EET: 7 _1 = 100 fs, J = 
50 cm' 1 , A = 100 cm -1 , T = 300 K , a regime in which the estimates in Table I do not 
apply. The initial excitation is assumed to be on site one. The Markovian approximation 
is seen to be very good for A = 1 cm -1 , fair for A = 2 cm -1 and invalid for reorganization 
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A = 1.0 cm" 1 A = 2.0 cm" 1 A = 10.0 cm" 1 




FIG. 2: Time evolution of population on site 1 [pn(i)] and the coherences [pi2(i)] [blue (solid) 
curve is the non-Mar kovian solution and red (dotted) curve is the Markovian approximation], for 
various values of A (in cm -1 ). Other parameter are: A = 100cm -1 , J = 50cm -1 , 7 = 10 13 sec -1 . 
The breakdown of the Markov approximation at A = 10 cm -1 is clearly visible. The time t' is 
dimensionless, with 10 units on this scale being equivalent to one ps. 

energies A > 10 cm -1 . 

To explore regimes of validity of the Markovian approximation for other values of the 
physical constants, we present sample results in Figs. [3] and SI obtained by varying (J, A) 
and (7 _1 , A), keeping A = 100 cm -1 . The results show for these cases that the Markovian 
approximation is poor for large A and "small' J and for large 7 _1 and small A. Other 
parameter values can be readily examined computationally using this approach. 
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J = 



100 cm" 1 , A = 2 cm" 1 



J = 60 cm 1 , A = 2 cm 1 



J = 5 cm 1 , A = 2 cm 1 




10 20 30 40 50 



100 cm" 1 , A = 10 cm" 1 





100 200 300 400 500 



FIG. 3: Time evolution of population on site 1: blue (solid) curve is the non-Markovian solution 
and red (dotted) curve is the Markovian approximation) for various values of the reorganization 
energy A and inter-site coupling J. The level separation A = 100 cm -1 and 7 _1 = 100 fs. It is 
clear that Markovian approximation is poor for large A and small J. Time t' is the dimensionless 
time, 10 units on this scale are equivalent to one ps. 

VI. REGIME OF VALIDITY OF THE SECOND-ORDER APPROXIMATION 



We have considered the master equation up to the second order in system-bath interac- 
tion. The aim of this section is to determine the system-bath interaction energy range (as 
represented by the reorganization energy A) over which the second order master equation [or 
coupled system of equations (Eq. (20)] can be used. To do so we compare estimates of the 
second order and the fourth order terms. This can be done analytically for the parameter 
regime where J ~ A. To do so we note that up to the second order, with the master equa- 
tion written in dimensionless time form [Eq. (20)], the magnitude of the second order term 



is of the order of 



for the case J ~ A. This arises by noting that \Hsr\ 2 ~ A, and the 
integral f* dr'e T '-*' cos[^(t' - t')]^') + rj 2 sin[^(t' - r')]y 2 (r / ) ~ 1, for the standard 
set of parameters (7- 1 = 100 fs, J = 20 cm" 1 , A = 100 cm' 1 , T = 300 K). 

Similarly, we can estimate the parameter dependence of the fourth order term. To esti- 
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= 1000 femtosec, A=2 cm 1 y '=200 femtosec, A=2 cm 1 y ' = 10 femtosec, A=2 cm 1 




20 40 60 80 100 



20 40 60 80 100 



100 200 300 400 500 



=1000 femtosec, A=10 cm" 1 y '=200 femtosec, A=10 cm 1 y '=10 femtosec, A=10 cm 1 




100 200 300 400 500 



FIG. 4: Time evolution of population on site 1 [blue (solid) curve is the non-Mar kovian solution 
and red (dotted) curve is the Markovian approximation] for various values of the reorganization 
energy A and phonon relaxation time 7 _1 . The level separation A = 100 cnr 1 and J = 50 cm . 
It is clear that Markovian approximation is poor for large relaxation times 7 -1 and small A, and 
it better for small 7 _1 . t' is the scaled time, as in the figure above. 

mate this we recall the Nakajima-Zwanzig master equation (valid to all orders) 
dp\t) 



dt 



j drtr R [c^Sit^QC^R^ p T (r). 



(29) 



where the time evolution operator is 

S(t,r) = T-+eM-ifdT , Q£ I SR (T')] 



(30) 



Here, Q = I — V is the well know projection operator and is the system-bath Liouvillian 
(—^[Hg R , .]). The time ordering operator orders time dependent operators from left to 
right with decreasing time arguments, to take into account the non-commutation of operators 
at different times. 

The zeroth order approximation to the time evolution operator S (t, r) gives the second 
order quantum master equation, the first order approximation to the time evolution operator 
gives the third order contribution which vanishes as the bath average of odd bath operators 
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vanish (see Ref. j^j), and the second order approximation to S gives the fourth order 
contribution. In order to estimate the magnitude of the latter term we write the fourth 
order term A 4 from Nakajima-Zwanzig equation as 

A *=^l dT f dr * i PsR®, Q&'sRfo), QiHUn), QiHU-r), Re q p(r)}}}} } ■ 

(31) 

We start from the interior commutator (1 — V)[Hg R (r), R eq p(r)} and recall that VO = 
ReqtrnO and Hg R = Vfu\ (where the summation convection is used). The bath average 
of single bath operators vanish [tr R (uj(r)R eq ) = 0] so that Q times the interior commu- 
tator gives [V/ (r)ui(r) , R eq p(r)]. Similarly, writing Q — 1 — V for the second from the 
interior commutator, and simplifying, the bath trace operation gives us the two time bath 
correlation functions {ui{r)ui{ri)) R . Repeating the same operations for the remaining com- 
mutators, noting that the bath averaging for the odd bath operators vanish and using the 
Wick theorem (uiUjU k Ui) R = (uiUj) R {u k u{) R + (uiU k ) R (ujui) R + (u^ui) R (ujU h ) R , we obtain 
that the fourth order term includes the product of two time bath correlation functions i.e., 
(ui(r)ui(Ti)) R (uj{Ti)Uj{j2)) R . On converting the equation into dimensionless time form as 
described in Section IV, and using the high temperature approximation for the correlation 
function, we conclude that the order-of-magnitude of the fourth order term is; 

A < - vh- (32) 

The ratio R42 of the fourth-order term to the second-order term is therefore R42 = 2^ ; 
which is of the order of 0.074 for A = 1 cm -1 , and ~ 0.74 for A = 10 cm" 1 . This suggests that 
the second order approximation for the master equation is good for A ~ 1 for the standard 
set of parameters (7^ = 100 fs, J = 20 cm -1 , A = 100 cm -1 , T = 300 K). However, 
for large A ~ 10 the fourth order term cannot be neglected. Interestingly, the domain of 
applicability of the second-order approximation in this J ~ A regime has a dependence on 
the same collection of parameters, (A / ' / y 2 /3h 2 small), as does the Markov approximation in 
the J < A. 
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Oscillator 2 



a 




Oscillator 1 



FIG. 5: Two harmonic oscillators excited by an ultra-short laser pulse. 
VII. INITIAL STATE PREPARATION BY AN ULTRA-SHORT LASER PULSE 
A. Formulation 

To understand the effect of initial coherences on the subsequent quantum dynamics, we 
consider a model of two 1-D harmonic oscillators separated by a distance a that are excited 
by an ultra short laser pulse (Fig. [5]). The results of this excitation are used below as 
sample coherent initial conditions for dimer propagation. Note that we restrict attention, as 
do most treatments of excitation of light-harvesting systems, to the "one-exciton manifold" , 
i.e. single excitations on each site. This model is useful to examine the dynamics, but should 
be augmented by excitation of states with bi-excitons of examining issues like entanglement, 
where the contributions from higher exciton states, no matter how small in magnitude, affect 
the entanglement measure. 

In the system (Fig. [5]) the wave vector of the laser pulse along the propagation direction 
makes an angle 9 with the line perpendicular that joining the oscillators. The laser field 
is treated semiclassically with the field sufficiently weak to allow first order perturbation 
theory for the light-oscillator interaction (22]. The total Hamiltonian in the coordinate 
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representation is then 



Ht = H svs + H, 



H S y s 



sys T 11 int 

h 2 d 2 l, 2 h 2 d 2 i 2 

2 + o^l - + o^|/ 2 



2m ch/^ 2 

ieh 



2m dy 2 ' 2 ' 



int 



2mc 



A.V 



A.V = A 1 (y 1 ,e,t)Lj— + A 2 (y 2 ,e,a,t)e3— , 



(33) 



with the vector potential A(y, 9, t) = ^ k e^(A k e^ kysine -^+Ale-^ kys ' ne -^). For a coherent 
laser pulse, e k = e. If A*{— ^) = ^(7), where c is the speed of light, then changing the 
summation to integration assuming continuous distribution of modes, we have 



/oo 
duA{u/c) exp[iu(y sm(6)/c — t)]. 
-00 



(34) 



For a Gaussian pulse A(u/c) = A exp[— i](ou — co^) 2 ], where uq is the central pulse fre- 
quency and i] defines the pulse width. The vector potentials take the form 

2" 



= A o\/- exp 



A 2 (yi,9,a,t) = A J-exp 



tu 



yi sin 9 



-t 



exp 



1 / yi sin 9 



4rj 



-t 



y 2 sm.9 a . 

— sin(0) - i 

c c 



exp 



1 ^2 sin a . 

sin(0) — i 

c c 



4r/ 



(35) 



The laser frequency is assumed tuned so as to excite the first excited state, with both 
oscillators initially (at t = —00) in their ground states. 
The eigensystems of oscillators 1 and 2 are 



(1) 




~~ 


vri/4 


(2) _ 




~~ 


vri/4 



1 

■— ( 
2 

1 

■-( 

2 



4" 

,( 2 ) 



V^2 3/2 r 1 2 21 



^ia 3 2 /2 y 2 exp[-^o£j/f](36) 



7T 



1/4' 



with o; C j = ^ki/m and a 4 = mki/h 2 , % = 1, 2. The total wavefunction of the system is 

${yi,y 2 ,t) = ^2a mn (t)u^ } (yi)v$ (y 2 ) exp[-^-(^ 1} + E$)t], a mn (t = -00) = 5 n0 5 m0 . 

(37) 
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Standard first-order perturbation theory gives the coefficients as 

a nm (t) = e o C ° S 6 I dt' [ dxi [ dx 2 u£ ) (x 1 )u^ ) (x 2 ) 
2mc J_ 00 J_ 00 

(d d \ 

A 1 (x 1 ,e,t')(—) + A 2 (x 2 ,e,a,t')(—)j u^^u^^eM-^nmt'} (38) 

Using the dipole approximation, the spatial integrals for both aoi and aio can be done 
exactly. The remaining time integral is treated as follows: the laser pulse is assumed to be 
ultrashort compared to the subsequent quantum dynamics. For times much greater than 
t/y/rj, the exponential e - ^ 2 / 4 ^ in time integration will be small, and the upper limit of the 
time integration can be extended to +00. The integration can then be performed exactly, 
giving 

a w = -Cv/^cosfle-^ 1 ^ 10 ) 2 , 

ooi = -C^e-^^) 2 . (39) 
yp 

Here \i = a 1 /a 2 = (h/k 2 ) 1/4 , ( = irA Q ey/aia 2 /\/2mc, u 01 = u 01 /uo, v w = oj w /u , with 
u mn = (E™ + - E$ - EP)/h. 

The first excited states of the both oscillators (E^\ E^) constitute our relevant system 
(they are separated by about 100 cm -1 in typical systems like FMO), and quantum dynamics 
takes place between them. The superposition of the excited states is written as 

|V) =ai |10)+ooi|01), (40) 

where 1 10} indicates that the first oscillator is excited and the second is in the ground state. 
The density matrix at the initial time is 

po = |VX^| = |ai | 2 |10)(10| + |a i| 2 |01)(01| + a 10 <|10>(01| + a 01 a* w \01) (10|. (41) 

This initial density matrix corresponds to a particular orientation angle 9. For excitation 
of an ensemble we average over theta, 

1 f 2n 

(po)e = ^J q PodO, (42) 

with the normalization 

P11 + P22 = (Ko| 2 )(? + (Ki| 2 )e = 1. (43) 
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k 2 



FIG. 6: (a) Population of angle averaged energy state |10)(10|, and (b) Coherence px2 as a function 
of fi = j^. Solid line, laser frequency ojq = 10 14 Hz, dashed line 5 x 10 14 Hz (Orange), dotted line 
10 15 Hz. Note that Imp 12 = 0. 



Thus, 



with 



(|aio| 2 )e 
(|ai | 2 )e + (|a i| 2 )e 

(|aoi| 2 )e 
(kio| 2 )e + (|aoi| 2 )e 
* _ (QioQpi)e 
9X2 ~ P21 ~ (\a w \ 2 )e + (\a i\ 2 )e' 



(ho|\ = icVe- 2 ^ {1+ ^ )2 , (|aoi| 2 ), = ^-e^ 1 ^ 



( a 10 a 0l)( 



_^2j i /^01^ e -^[(l +no )2+(l +i/01 )2 ] _ 



(44) 



(45) 



aujQ\ ' ' c 

Here Ji(.) is the Bessel function of first kind and of order 1. This constitutes the initial 
density matrix for the relevant system. 



B. Numerical Results 

Let 7] = 1CT 4 ps 2 , the separation between the two oscillators a = 2 nm, Co> i = —^ C 2 — 
— 10 14 Hz, and ojiq = v //Iwoi- Figure |6] shows the density matrix elements [i.e. the initial 
state given by Eqs. (144j) and (145p ] as a function of the ratio of oscillator for constants kx/k^ 
for various values of the excitation laser frequency. 

This initial state can then be used as a model for the initial conditions for the numerical 
solution of the non-Markovian equations [Eq. ( !T5|) ] for the dimer. Figure [7] (in the site 
representation) compares the time evolution of pxx(t) for an initial "no-coherence" state 
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ki/jfe 2 = 1.0 kjk 2 = 0.5 ^2 = 0.1 



Pn Pn Pn 




FIG. 7: The effect of initial coherence on the of relaxation dynamics. Long time behavior for 
various values of ki/hz- Here A = 1 cm -1 and ujq = 10 14 . The blue (solid) curve is for "no- 
coherence" initial condition , and red (dotted) curve is for the "initial coherence" condition, i.e., 
initial value of the density matrix elements for the model photo-excitation studied above. The 
second and third row shows the dynamics of the off-diagonal elements. 

(only populations) and the model generated state with "initial coherence" , for various values 
of ki/k 2 . 

One may identify two time scales associated with the electronic energy transfer, the time 
scale over which the site occupation pn becomes relatively constant, and the rate at which 
this occurs. From the plots of pn in Fig. [7J it appears that the presence of initial coherence 
( at t — 0) effects both of these time scales, but has little effect on the overall damping-out 
of the coherence, i.e. the overall decay of oscillations in both the real and imaginary parts of 
pi 2- By contrast, the fall-off rate for the decay of pn is far faster for the case with initially 



21 



no-coherence than it is for the case where there initially is coherence. In cases other than 
k\/k2 = 1.0 the time at which the system reaches the equilibrium value of 1/2 seems similar 
in both the cases where there is coherence initially and where there is not. 

VIII. CONCLUSION 

A straightforward approach to solving the second order Born master equation, with and 
without the Markov approximation, has been introduced. In addition to obtaining numer- 
ical results showing the range of validity of these approximations, a number of analytical 
estimates, shown in Table I, of parameter ranges over which these approximations can be 
used has been obtained. For the case of the traditional dimer model for electronic energy 
transfer in photosynthesis, surprisingly small reorganization energies (a few cm -1 ) are re- 
quired for the validity of the Markovian approximation. In addition, we note that for dimer 
coupling strengths on the order of the energy difference between site energies, higher order 
terms than second order in the system-bath coupling are required if AX/(h 2 /3'y 2 ) « 1 is not 
satisfied, where A is the reorganization energy, and 7 defines the exponential falloff rate of 
the bath correlation function. Once again, the limitation to small reorganization energies, 
not well appreciated in the past, is made explicit. 

We have also provided an example of the role of initial coherences in the subsequent 
evolution of the dimer dynamics for typical parameters associated with model photosynthetic 
light harvesting systems. 
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